{
 "cells": [
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "# 遥感图像分类 - 随机森林方法",
   "id": "83cb767a4646e34f"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 1. 环境准备与数据加载",
   "id": "42b77f320143c84e"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import numpy as np\n",
    "import rasterio\n",
    "from rasterio.plot import show\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import confusion_matrix, accuracy_score, cohen_kappa_score\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "from osgeo import gdal\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "\n",
    "# 设置中文显示\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']\n",
    "plt.rcParams['axes.unicode_minus'] = False"
   ],
   "id": "a6db6071b023431f",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 2. 参数设置",
   "id": "1dfe666ec4696d9b"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 输入输出路径设置\n",
    "input_image = \"../data/output_preprocessed/preprocessed_landsat8_2021.tif\"  # 预处理后的影像\n",
    "training_samples = \"../data/training_samples/auto_samples/auto_samples_xinDu_2021.shp\"  # 训练样本\n",
    "output_dir = \"../data/output_classification\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "\n",
    "raw_result = os.path.join(output_dir, \"landcover_xinDu_2021_raw.tif\")  # 原始分类结果\n",
    "postprocessed_result = os.path.join(output_dir, \"landcover_xinDu_2021_post.tif\")  # 后处理结果\n",
    "\n",
    "# 类别定义 (需要与样本shapefile中的class字段一致)\n",
    "CLASS_NAMES = {\n",
    "    1: \"林地\",\n",
    "    2: \"草地\",\n",
    "    3: \"城镇区域\",\n",
    "    4: \"水域\",\n",
    "    5: \"灌木林地\",\n",
    "    6: \"农田\",\n",
    "    7: \"裸地\",\n",
    "    8: \"其他\"\n",
    "}\n",
    "\n",
    "# 类别颜色 (RGB)\n",
    "CLASS_COLORS = {\n",
    "    1: (0, 100, 0),  # 林地 - 深绿\n",
    "    2: (154, 205, 50),  # 草地 - 黄绿\n",
    "    3: (178, 34, 34),  # 城镇 - 砖红\n",
    "    4: (65, 105, 225),  # 水域 - 蓝色\n",
    "    5: (34, 139, 34),  # 灌木 - 绿色\n",
    "    6: (255, 215, 0),  # 农田 - 金色\n",
    "    7: (139, 137, 137),  # 裸地 - 灰色\n",
    "    8: (255, 255, 255)  # 其他 - 白色\n",
    "}"
   ],
   "id": "4f12fd8b8f590bd5",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# CLASS_NAMES\n",
    "CLASS_COLORS"
   ],
   "id": "78daf89a23cc9a90",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 3. 加载训练样本",
   "id": "45cbfbe9c7c30b49"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def load_training_samples(sample_shp, raster_path):\n",
    "    \"\"\"\n",
    "    加载训练样本数据\n",
    "    :param sample_shp: 样本点shapefile路径\n",
    "    :param raster_path: 对应的栅格图像路径\n",
    "    :return: 特征矩阵X和标签y\n",
    "    \"\"\"\n",
    "    print(\"正在加载训练样本...\")\n",
    "    # 读取样本点\n",
    "    samples = gpd.read_file(sample_shp)\n",
    "\n",
    "    # 读取栅格数据\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        # 提取样本点对应的像素值\n",
    "        X = []\n",
    "        y_labels = []  # 改名为y_labels避免冲突\n",
    "\n",
    "        # 使用tqdm显示进度条\n",
    "        for idx, row in tqdm(samples.iterrows(), total=len(samples), desc=\"提取样本点数据\"):\n",
    "            # 获取样本点坐标\n",
    "            x, y_coord = row.geometry.x, row.geometry.y\n",
    "            # 将地理坐标转换为图像坐标\n",
    "            row_idx, col_idx = src.index(x, y_coord)\n",
    "            # 读取该位置所有波段的值\n",
    "            pixel_values = src.read()[:, row_idx, col_idx]\n",
    "            X.append(pixel_values)\n",
    "            y_labels.append(row['class'])  # 确保shapefile中有'class'字段\n",
    "\n",
    "    return np.array(X), np.array(y_labels)\n",
    "\n",
    "\n",
    "# 示例调用\n",
    "training_samples = \"../data/training_samples/auto_samples/auto_samples_xinDu_2021.shp\"\n",
    "input_image = \"../data/output_preprocessed/preprocessed_landsat8_2021.tif\"\n",
    "X, y = load_training_samples(training_samples, input_image)\n",
    "\n"
   ],
   "id": "6bb850a2b054789b",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 可视化样本统计信息",
   "id": "dd119cc47fd4209f"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 显示样本统计信息\n",
    "print(\"\\n样本统计信息:\")\n",
    "\n",
    "sample_counts = pd.Series(y).value_counts().sort_index()\n",
    "# 类别名称映射，可以这样显示（假设CLASS_NAMES已定义）\n",
    "sample_counts.index = [CLASS_NAMES[c] for c in sample_counts.index]\n",
    "print(sample_counts)"
   ],
   "id": "8c42be92b32df2cc",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 1. 加载数据\n",
    "X, y = load_training_samples(training_samples, input_image)\n",
    "\n",
    "# 2. 训练模型\n",
    "model = RandomForestClassifier(n_estimators=100, random_state=42)\n",
    "model.fit(X, y)\n",
    "\n",
    "# 3. 预测整个影像\n",
    "with rasterio.open(input_image) as src:\n",
    "    img = src.read()\n",
    "    height, width = src.height, src.width\n",
    "    pixels = img.reshape(img.shape[0], -1).T\n",
    "\n",
    "# 分块预测\n",
    "y_pred_full = np.zeros(height * width, dtype=np.uint8)\n",
    "chunk_size = 10000\n",
    "for i in tqdm(range(0, len(pixels), chunk_size)):\n",
    "    y_pred_full[i:i + chunk_size] = model.predict(pixels[i:i + chunk_size])\n",
    "y_pred_2d = y_pred_full.reshape(height, width)\n",
    "\n",
    "# 4. 生成彩色结果\n",
    "colored_result = np.zeros((height, width, 3), dtype=np.uint8)\n",
    "for class_id, color in CLASS_COLORS.items():\n",
    "    colored_result[y_pred_2d == class_id] = color\n",
    "\n",
    "# 5. 可视化\n",
    "with rasterio.open(input_image) as src:\n",
    "    rgb = src.read([4, 3, 2])  # 假彩色合成\n",
    "    rgb = np.transpose(rgb, (1, 2, 0))\n",
    "    rgb = (rgb - rgb.min()) / (rgb.max() - rgb.min())\n",
    "\n",
    "patches = [plt.Rectangle((0, 0), 1, 1, color=np.array(color) / 255)\n",
    "           for color in CLASS_COLORS.values()]\n",
    "\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.imshow(rgb)\n",
    "plt.imshow(colored_result, alpha=0.5)\n",
    "plt.title(\"分类结果叠加原始影像\", fontsize=16)\n",
    "plt.axis('off')\n",
    "plt.legend(patches, CLASS_NAMES.values(), bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(output_dir, \"overlay_map.png\"), dpi=300, bbox_inches='tight')\n",
    "plt.show()"
   ],
   "id": "2078be898aabb468",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 4. 训练随机森林分类器\n",
    "- stratify=y：按类别分层抽样，保持训练集和测试集中各类别比例与原数据一致（避免某些类别样本过少）\n",
    "- 随机森林是一种集成学习算法，适合处理高维遥感数据，能自动选择重要特征并减少过拟合。参数详解：\n",
    "    - `n_estimators=100`：使用100棵决策树（树越多效果通常越好，但计算量增大）。\n",
    "    - `random_state=42`：固定随机种子，确保模型可复现。\n",
    "    - `n_jobs=-1`：使用所有CPU核心并行加速训练。\n",
    "    - `class_weight='balanced'`：自动平衡类别权重，解决类别样本不均衡问题（如城镇样本远少于林地时特别重要）。\n",
    "    - `oob_score=True`：计算袋外误差（Out-of-Bag Score），用于快速评估模型泛化能力（无需单独测试集）。\n",
    "- 为什么选择随机森林？\n",
    "    1. **抗噪声能力强**：适合处理遥感数据中的混合像元问题。\n",
    "    2. **自动特征选择**：通过`feature_importances_`可知哪些波段对分类最重要（如近红外比红光更重要）。\n",
    "    3. **并行计算**：`n_jobs=-1`能显著加速大规模数据训练。\n",
    "\n"
   ],
   "id": "9ceab01a9df422bf"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 划分训练集和测试集\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X, y, test_size=0.3, random_state=42, stratify=y)\n",
    "\n",
    "# 训练随机森林分类器\n",
    "print(\"\\n训练随机森林分类器...\")\n",
    "rf = RandomForestClassifier(n_estimators=100,\n",
    "                            random_state=42,\n",
    "                            n_jobs=-1,\n",
    "                            class_weight='balanced',\n",
    "                            oob_score=True)\n",
    "rf.fit(X_train, y_train)"
   ],
   "id": "49528543e3b316cf",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "### 输出随机森林分类器的关键评估指标\n",
    "\n",
    "1. `rf.oob_score_`—— 袋外分数 (Out-of-Bag Score)\n",
    "   - **作用**：评估模型泛化能力的指标（类似测试集准确率）。\n",
    "   - **原理**：\n",
    "     - 随机森林每棵树训练时，随机抽取约63%的样本，剩余37%未参与训练的样本称为\"袋外样本\"（OOB）。\n",
    "     - 用这些袋外样本评估每棵树的预测效果，最终汇总所有树的OOB预测结果计算整体准确率。\n",
    "   - **解读**：\n",
    "     - 值范围：`0~1`，越接近1表示模型越好。\n",
    "     - 例如 `OOB分数: 0.85`表示模型对未知数据的预测准确率约85%。\n",
    "   - **优点**：无需额外划分测试集即可快速评估模型。\n",
    "2. `rf.feature_importances_`—— 特征重要性\n",
    "   * **作用**：分析各波段（特征）对分类结果的贡献度。\n",
    "   * **原理**：\n",
    "     - 通过计算每个波段在决策树分裂时的信息增益（Gini不纯度减少或熵减少），统计所有树中的平均贡献。\n",
    "     - 结果是一个数组，长度等于输入特征数（如Landsat 8的7个波段对应7个值）。"
   ],
   "id": "6e926c9a2dbc3b97"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 打印模型信息\n",
    "print(f\"OOB分数: {rf.oob_score_:.3f}\")  # 袋外估计分数\n",
    "print(f\"特征重要性: {rf.feature_importances_}\")      # 查看各波段对分类的贡献"
   ],
   "id": "aeab857772e59fa2",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 5. 模型评估\n",
    "\n",
    "### 遥感影像分类结果的综合评估工具\n",
    "> 用于量化分类模型的精度并生成可视化报告。\n",
    "\n",
    "\n",
    "\n",
    "#### 1. 核心功能\n",
    "\n",
    "- **输入**：\n",
    "  - `y_true`：真实的地物类别标签（如测试集真实值）\n",
    "  - `y_pred`：模型预测的类别标签\n",
    "  - `class_names`：类别名称列表（如 `[\"林地\", \"草地\", ...]`）\n",
    "- **输出**：\n",
    "  - 总体精度、Kappa系数、各类别精度表格\n",
    "  - 归一化混淆矩阵热力图\n",
    "\n",
    "#### 2. 关键评估指标\n",
    "\n",
    "##### (1) 总体精度 (Overall Accuracy)\n",
    "\n",
    "```\n",
    "overall_accuracy = accuracy_score(y_true, y_pred)\n",
    "```\n",
    "\n",
    "- **含义**：所有正确分类的样本占总样本的比例。\n",
    "- **示例**：`0.85`表示85%的像素被正确分类。\n",
    "- **局限性**：对样本不均衡的数据可能误导（如90%是林地时，全预测为林地也可达90%准确率）。\n",
    "\n",
    "##### (2) Kappa系数 (Cohen's Kappa)\n",
    "\n",
    "```\n",
    "kappa = cohen_kappa_score(y_true, y_pred)\n",
    "```\n",
    "\n",
    "- **含义**：剔除随机分类影响后的分类一致性，范围 `[-1, 1]`。\n",
    "- **解读**：\n",
    "  - `<0.4`：差\n",
    "  - `0.4-0.75`：中到好\n",
    "  - `>0.75`：优秀\n",
    "- **用途**：比总体精度更能反映模型真实能力。\n",
    "\n",
    "##### (3) 生产者精度 & 用户精度\n",
    "\n",
    "```\n",
    "cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n",
    "class_report = pd.DataFrame({\n",
    "    '生产者精度': cm.diagonal(),  # 每类被正确分类的比例\n",
    "    '用户精度': cm.T.diagonal()   # 预测为某类的样本中实际正确的比例\n",
    "})\n",
    "```\n",
    "\n",
    "- **生产者精度（Producer's Accuracy）**：\n",
    "  - 真实为“林地”的像素有多少被正确分类。\n",
    "  - **公式**：`正确分类的林地像素 / 真实林地总像素`\n",
    "- **用户精度（User's Accuracy）**：\n",
    "  - 模型预测为“林地”的像素中有多少真是林地。\n",
    "  - **公式**：`正确分类的林地像素 / 模型预测的所有林地像素`\n",
    "\n",
    "#### 3. 混淆矩阵可视化\n",
    "\n",
    "```\n",
    "plt.imshow(cm, cmap=plt.cm.Blues)\n",
    "```\n",
    "\n",
    "- **作用**：直观展示各类别的错分情况。\n",
    "- **示例解读**：\n",
    "  - 行：真实类别\n",
    "  - 列：预测类别\n",
    "  - 对角线越亮（值接近1）表示分类越好\n",
    "  - 非对角线亮斑表示常见混淆（如“裸地”被误分为“城镇”）"
   ],
   "id": "d5ba00f25023142b"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def evaluate_classification(y_true, y_pred, class_names):\n",
    "    \"\"\"\n",
    "    评估分类精度\n",
    "    :param y_true: 真实标签\n",
    "    :param y_pred: 预测标签\n",
    "    :param class_names: 类别名称列表\n",
    "    \"\"\"\n",
    "    # 计算混淆矩阵\n",
    "    cm = confusion_matrix(y_true, y_pred)\n",
    "\n",
    "    # 计算总体精度\n",
    "    overall_accuracy = accuracy_score(y_true, y_pred)\n",
    "\n",
    "    # 计算Kappa系数\n",
    "    kappa = cohen_kappa_score(y_true, y_pred)\n",
    "\n",
    "    # 打印结果\n",
    "    print(\"\\n分类精度评估:\")\n",
    "    print(f\"总体精度: {overall_accuracy:.2%}\")\n",
    "    print(f\"Kappa系数: {kappa:.4f}\")\n",
    "\n",
    "    # 计算并打印各类别生产者精度和用户精度\n",
    "    print(\"\\n各类别精度:\")\n",
    "    cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]  # 归一化\n",
    "    class_report = pd.DataFrame({\n",
    "        '生产者精度': cm.diagonal(),\n",
    "        '用户精度': cm.T.diagonal()\n",
    "    }, index=class_names)\n",
    "    print(class_report.round(4))\n",
    "\n",
    "    # 可视化混淆矩阵\n",
    "    plt.figure(figsize=(10, 8))\n",
    "    plt.imshow(cm, interpolation='nearest', cmap=plt.cm.Blues)\n",
    "    plt.title(\"混淆矩阵 (归一化)\")\n",
    "    plt.colorbar()\n",
    "    tick_marks = np.arange(len(class_names))\n",
    "    plt.xticks(tick_marks, class_names, rotation=45)\n",
    "    plt.yticks(tick_marks, class_names)\n",
    "    plt.ylabel('真实类别')\n",
    "    plt.xlabel('预测类别')\n",
    "\n",
    "    # 在矩阵中显示数值\n",
    "    thresh = cm.max() / 2.\n",
    "    for i in range(cm.shape[0]):\n",
    "        for j in range(cm.shape[1]):\n",
    "            plt.text(j, i, f\"{cm[i, j]:.2f}\",\n",
    "                    horizontalalignment=\"center\",\n",
    "                    color=\"white\" if cm[i, j] > thresh else \"black\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "# 在测试集上评估\n",
    "print(\"在测试集上评估模型...\")\n",
    "y_pred = rf.predict(X_test)\n",
    "evaluate_classification(y_test, y_pred, list(CLASS_NAMES.values()))"
   ],
   "id": "cb75b4714c540886",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": [
    "## 6. 影像分类\n",
    "\n",
    "这段代码是一个**用于对整幅遥感影像进行分类的Python函数**，它使用训练好的机器学习模型（如随机森林）对输入的多波段遥感影像进行逐像素分类，并将分类结果保存为新的栅格文件。==通过这段代码，高效地将机器学习模型应用于大面积遥感影像分类，并生成可直接在GIS软件（如QGIS）中使用的结果。==\n",
    "\n",
    "### 1、函数功能\n",
    "\n",
    "- **输入**：\n",
    "  - `raster_path`：待分类的遥感影像路径（如Landsat 8的TIFF文件）\n",
    "  - `model`：训练好的分类模型（如`RandomForestClassifier`）\n",
    "  - `output_path`：分类结果输出路径\n",
    "- **输出**：\n",
    "  - 保存分类结果为单波段TIFF文件，每个像素值代表预测的类别标签（如1=林地，2=草地等）\n",
    "\n",
    "### 2、为什么分块？\n",
    "\n",
    "- 大影像（如10000×10000像素）直接读取会耗尽内存，分块处理（如512×512）可降低内存需求。\n",
    "- `classification`矩阵用于存储最终的分类结果。\n",
    "- **关键步骤**：\n",
    "  1. **窗口读取**：`src.read(window=window)`只读取当前块的像素数据。\n",
    "  2. **数据重塑**：将 `(bands, height, width)`转为 `(n_pixels, n_bands)`，符合sklearn输入格式。\n",
    "  3. **模型预测**：对当前块的所有像素进行分类。\n",
    "  4. **结果填充**：将预测结果放回 `classification`矩阵的对应位置。\n",
    "\n",
    "### 3. 使用场景\n",
    "\n",
    "1. **土地利用分类**：将影像分为林地、水域、城镇等类别。\n",
    "2. **变化检测**：对不同年份的影像分类后比较变化区域。\n",
    "3. **生态研究**：分析植被覆盖动态。"
   ],
   "id": "4b06cee0de9544fc"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def classify_image(raster_path, model, output_path):\n",
    "    \"\"\"\n",
    "    对整个图像进行分类\n",
    "    :param raster_path: 输入栅格路径\n",
    "    :param model: 训练好的分类模型\n",
    "    :param output_path: 分类结果输出路径\n",
    "    \"\"\"\n",
    "    print(\"\\n开始影像分类...\")\n",
    "    # 读取栅格数据\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        # 读取所有波段数据和元数据\n",
    "        profile = src.profile\n",
    "        rows, cols = src.shape\n",
    "        n_bands = src.count\n",
    "\n",
    "        # 分块处理大影像 (避免内存不足)\n",
    "        block_size = 512  # 分块大小\n",
    "        classification = np.zeros((rows, cols), dtype=np.uint8)\n",
    "\n",
    "        # 计算分块数量\n",
    "        n_blocks_x = int(np.ceil(cols / block_size))\n",
    "        n_blocks_y = int(np.ceil(rows / block_size))\n",
    "\n",
    "        # 分块处理\n",
    "        for i in tqdm(range(n_blocks_y), desc=\"处理行块\"):\n",
    "            for j in range(n_blocks_x):\n",
    "                # 计算当前块的起始和结束位置\n",
    "                y_start = i * block_size\n",
    "                y_end = min((i + 1) * block_size, rows)\n",
    "                x_start = j * block_size\n",
    "                x_end = min((j + 1) * block_size, cols)\n",
    "\n",
    "                # 读取当前块的数据\n",
    "                window = ((y_start, y_end), (x_start, x_end))\n",
    "                block_data = src.read(window=window)\n",
    "\n",
    "                # 重塑数据为(像素数, 波段数)\n",
    "                block_shape = block_data.shape\n",
    "                reshaped_data = block_data.reshape(block_shape[0], -1).T\n",
    "\n",
    "                # 预测当前块\n",
    "                block_pred = model.predict(reshaped_data)\n",
    "\n",
    "                # 将预测结果放回正确位置\n",
    "                classification[y_start:y_end, x_start:x_end] = block_pred.reshape(\n",
    "                    y_end - y_start, x_end - x_start)\n",
    "\n",
    "        # 更新元数据\n",
    "        profile.update(\n",
    "            dtype=rasterio.uint8,\n",
    "            count=1,\n",
    "            nodata=0,\n",
    "            compress='lzw'  # 使用LZW压缩减小文件大小\n",
    "        )\n",
    "\n",
    "        # 保存结果\n",
    "        with rasterio.open(output_path, 'w', **profile) as dst:\n",
    "            dst.write(classification, 1)\n",
    "\n",
    "    print(f\"分类结果已保存至: {output_path}\")\n",
    "\n",
    "# 执行分类\n",
    "classify_image(input_image, rf, raw_result)"
   ],
   "id": "3d849d754927388d",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 可视化结果",
   "id": "7f9270d621dc321c"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 读取分类结果\n",
    "with rasterio.open(raw_result) as src:\n",
    "    classified = src.read(1)\n",
    "\n",
    "# 生成彩色图像\n",
    "colors = np.array(list(CLASS_COLORS.values())) / 255\n",
    "colored_map = colors[classified - 1]  # 类别1对应colors[0]\n",
    "\n",
    "# 显示\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.imshow(colored_map)\n",
    "plt.legend(handles=[\n",
    "    plt.Rectangle((0,0),1,1, color=color) for color in colors\n",
    "], labels=CLASS_NAMES.values())\n",
    "plt.title(\"分类结果\", fontsize=16)\n",
    "plt.axis('off')\n",
    "plt.show()"
   ],
   "id": "289c9f4b55b23e44",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 7. 分类后处理",
   "id": "3cc6161b3109d6a3"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def postprocess_classification(input_path, output_path, min_size=4):\n",
    "    \"\"\"\n",
    "    分类后处理 - 去除小图斑\n",
    "    :param input_path: 分类结果路径\n",
    "    :param output_path: 后处理结果路径\n",
    "    :param min_size: 最小图斑大小(像素数)\n",
    "    \"\"\"\n",
    "    print(\"\\n执行分类后处理...\")\n",
    "    # 使用GDAL进行聚类处理\n",
    "    src_ds = gdal.Open(input_path)\n",
    "    band = src_ds.GetRasterBand(1)\n",
    "\n",
    "    # 创建输出文件\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    dst_ds = driver.Create(output_path,\n",
    "                          src_ds.RasterXSize,\n",
    "                          src_ds.RasterYSize,\n",
    "                          1,\n",
    "                          band.DataType)\n",
    "    dst_ds.SetGeoTransform(src_ds.GetGeoTransform())\n",
    "    dst_ds.SetProjection(src_ds.GetProjection())\n",
    "\n",
    "    # 执行聚类\n",
    "    gdal.SieveFilter(band, None, dst_ds.GetRasterBand(1),\n",
    "                    min_size, 4)\n",
    "\n",
    "    # 关闭数据集\n",
    "    src_ds = None\n",
    "    dst_ds = None\n",
    "\n",
    "    print(f\"后处理结果已保存至: {output_path}\")\n",
    "\n",
    "# 执行后处理\n",
    "postprocess_classification(raw_result, postprocessed_result, min_size=10)"
   ],
   "id": "d67e5f629e3232a9",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 可视化结果",
   "id": "3e01cb3e24d09155"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 读取分类结果\n",
    "with rasterio.open(postprocessed_result) as src:\n",
    "    classified = src.read(1)\n",
    "\n",
    "# 生成彩色图像\n",
    "colors = np.array(list(CLASS_COLORS.values())) / 255\n",
    "colored_map = colors[classified - 1]  # 类别1对应colors[0]\n",
    "\n",
    "# 显示\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.imshow(colored_map)\n",
    "plt.legend(handles=[\n",
    "    plt.Rectangle((0,0),1,1, color=color) for color in colors\n",
    "], labels=CLASS_NAMES.values())\n",
    "plt.title(\"分类结果\", fontsize=16)\n",
    "plt.axis('off')\n",
    "plt.show()"
   ],
   "id": "976cecf63338e211",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "### 对比",
   "id": "122ef09b638d0d21"
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "#### 1.并排对比法",
   "id": "32160a2d72f1e1b6"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import rasterio\n",
    "\n",
    "# 读取两个结果\n",
    "with rasterio.open(raw_result) as src:\n",
    "    raw_classified = src.read(1)\n",
    "    profile = src.profile  # 保存元数据用于后续分析\n",
    "\n",
    "with rasterio.open(postprocessed_result) as src:\n",
    "    post_classified = src.read(1)\n",
    "\n",
    "# 生成颜色映射\n",
    "colors = np.array(list(CLASS_COLORS.values())) / 255\n",
    "\n",
    "# 创建对比图\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))\n",
    "\n",
    "# 原始结果\n",
    "ax1.imshow(colors[raw_classified - 1])\n",
    "ax1.set_title(\"原始分类结果\", fontsize=14)\n",
    "ax1.axis('off')\n",
    "\n",
    "# 后处理结果\n",
    "im = ax2.imshow(colors[post_classified - 1])\n",
    "ax2.set_title(\"后处理结果（去小斑块）\", fontsize=14)\n",
    "ax2.axis('off')\n",
    "\n",
    "# 添加公共图例\n",
    "fig.legend(\n",
    "    handles=[plt.Rectangle((0,0),1,1, color=color) for color in colors],\n",
    "    labels=CLASS_NAMES.values(),\n",
    "    loc='lower center',\n",
    "    ncol=len(CLASS_NAMES),\n",
    "    bbox_to_anchor=(0.5, -0.05)\n",
    ")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(output_dir, \"comparison.png\"), dpi=300, bbox_inches='tight')\n",
    "plt.show()"
   ],
   "id": "fff67746c2d24e35",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "#### 2. 差异高亮法（显示变化区域）",
   "id": "ddfa1c85536c65e2"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 计算差异区域 (0表示无变化，1表示被修改)\n",
    "diff = (raw_classified != post_classified).astype(np.uint8)\n",
    "\n",
    "plt.figure(figsize=(12, 10))\n",
    "plt.imshow(colors[post_classified - 1])  # 显示后处理结果为底图\n",
    "plt.imshow(diff, cmap='Reds', alpha=0.3)  # 红色半透明叠加差异区域\n",
    "plt.title(\"后处理修改区域（红色标记）\", fontsize=16)\n",
    "plt.axis('off')\n",
    "\n",
    "# 添加图例\n",
    "from matplotlib.lines import Line2D\n",
    "legend_elements = [\n",
    "    Line2D([0], [0], marker='s', color='w', label='被修改区域',\n",
    "           markerfacecolor='red', markersize=15)\n",
    "]\n",
    "plt.legend(handles=legend_elements, loc='lower right')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(output_dir, \"changes_highlight.png\"), dpi=300)\n",
    "plt.show()"
   ],
   "id": "56c308b34b3f5469",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "#### 3. 分屏交互对比（高级）",
   "id": "5f88ad8fcdeab6ad"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "from mpl_toolkits.axes_grid1 import ImageGrid\n",
    "\n",
    "fig = plt.figure(figsize=(20, 10))\n",
    "grid = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.3)\n",
    "\n",
    "grid[0].imshow(colors[raw_classified - 1])\n",
    "grid[0].set_title(\"原始分类结果\")\n",
    "grid[1].imshow(colors[post_classified - 1])\n",
    "grid[1].set_title(\"后处理结果\")\n",
    "\n",
    "# 同步缩放和移动\n",
    "def on_move(event):\n",
    "    if event.inaxes == grid[0]:\n",
    "        grid[1].set_xlim(grid[0].get_xlim())\n",
    "        grid[1].set_ylim(grid[0].get_ylim())\n",
    "    elif event.inaxes == grid[1]:\n",
    "        grid[0].set_xlim(grid[1].get_xlim())\n",
    "        grid[0].set_ylim(grid[1].get_ylim())\n",
    "\n",
    "fig.canvas.mpl_connect('motion_notify_event', on_move)\n",
    "plt.show()"
   ],
   "id": "500bac3d6461f3ae",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 8. 结果可视化",
   "id": "d2487e69eb611751"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "def visualize_classification(result_path, class_colors):\n",
    "    \"\"\"\n",
    "    可视化分类结果\n",
    "    :param result_path: 分类结果路径\n",
    "    :param class_colors: 类别颜色字典 {类别ID: RGB颜色}\n",
    "    \"\"\"\n",
    "    print(\"\\n可视化分类结果...\")\n",
    "    with rasterio.open(result_path) as src:\n",
    "        data = src.read(1)\n",
    "\n",
    "        # 创建彩色图像\n",
    "        rgb = np.zeros((data.shape[0], data.shape[1], 3), dtype=np.uint8)\n",
    "\n",
    "        for class_id, color in class_colors.items():\n",
    "            rgb[data == class_id] = color\n",
    "\n",
    "        # 显示结果\n",
    "        plt.figure(figsize=(12, 10))\n",
    "        plt.imshow(rgb)\n",
    "        plt.title(\"土地覆盖分类结果\")\n",
    "        plt.axis('off')\n",
    "\n",
    "        # 创建图例\n",
    "        patches = [plt.Rectangle((0,0),1,1, color=np.array(color)/255)\n",
    "                   for color in class_colors.values()]\n",
    "        plt.legend(patches, CLASS_NAMES.values(),\n",
    "                  bbox_to_anchor=(1.05, 1),\n",
    "                  loc='upper left',\n",
    "                  title=\"土地覆盖类别\")\n",
    "\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "# 可视化原始分类结果和后处理结果\n",
    "print(\"原始分类结果:\")\n",
    "visualize_classification(raw_result, CLASS_COLORS)\n",
    "\n",
    "print(\"\\n后处理结果:\")\n",
    "visualize_classification(postprocessed_result, CLASS_COLORS)"
   ],
   "id": "9d2a3d43839197a6",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "markdown",
   "source": "## 9. 保存分类报告",
   "id": "387f511f7332389f"
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": [
    "# 生成分类报告\n",
    "def generate_classification_report(model, output_path):\n",
    "    \"\"\"\n",
    "    生成分类报告\n",
    "    :param model: 训练好的模型\n",
    "    :param output_path: 报告保存路径\n",
    "    \"\"\"\n",
    "    # 获取特征重要性\n",
    "    feature_importances = model.feature_importances_\n",
    "\n",
    "    # 创建报告\n",
    "    report = {\n",
    "        \"总体精度\": accuracy_score(y_test, y_pred),\n",
    "        \"Kappa系数\": cohen_kappa_score(y_test, y_pred),\n",
    "        \"OOB分数\": model.oob_score_,\n",
    "        \"特征重要性\": feature_importances,\n",
    "        \"类别分布\": pd.Series(y).value_counts().sort_index().to_dict()\n",
    "    }\n",
    "\n",
    "    # 保存报告\n",
    "    with open(output_path, 'w', encoding='utf-8') as f:\n",
    "        f.write(\"=== 分类报告 ===\\n\\n\")\n",
    "        f.write(f\"总体精度: {report['总体精度']:.2%}\\n\")\n",
    "        f.write(f\"Kappa系数: {report['Kappa系数']:.4f}\\n\")\n",
    "        f.write(f\"OOB分数: {report['OOB分数']:.4f}\\n\\n\")\n",
    "\n",
    "        f.write(\"特征重要性:\\n\")\n",
    "        for i, imp in enumerate(report['特征重要性']):\n",
    "            f.write(f\"波段 {i+1}: {imp:.4f}\\n\")\n",
    "\n",
    "        f.write(\"\\n类别分布:\\n\")\n",
    "        for class_id, count in report['类别分布'].items():\n",
    "            f.write(f\"{CLASS_NAMES[class_id]}: {count} 样本\\n\")\n",
    "\n",
    "    print(f\"分类报告已保存至: {output_path}\")\n",
    "\n",
    "# 生成报告\n",
    "generate_classification_report(rf, os.path.join(output_dir, \"classification_report.txt\"))"
   ],
   "id": "edd011c747710724",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "d49be4c68fd86d79",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "b76013f7d1394858",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "92aeee3c4ed52f6",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "id": "8a405c00bc69cbb6",
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
